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Shape dependence of higher order correlations introduces complication in direct determination of these 
quantities. For this reason theoretical and observational progress has been restricted in calculating one 
point distribution functions and related moments. Methods based on factorial moments of two point count 
probability distribution (cumulant correlators) were recently shown to be efficient in subtracting discreteness 
effects and extracting useful information from galaxy catalogs. We use these cumulant correlators Qnm to 
study clustering in scale free simulations both in two and three dimensions. Using this method in the 
highly nonlinear regime we were able to separate hierarchal amplitudes Q^,R a and Rb associated with 
different tree graphs contributing to third and fourth order correlation functions. They were found to 
increase with power on large scales. Results based on factorial moments of one point cumulants show very 
good agreement with that on cumulant correlators. Comparison of simulation results with perturbation 
theory and extended perturbation theory were found to be in reasonable agreement. Comparisons were also 
made against predictions from hierarchal models for higher order correlations. We argue that finite volume 
corrections are very important for computation of cumulant correlators. 
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1. Introduction 

2 ' There is growing evidence which suggest that the large scale structure in the universe was formed by the 
\ gravitational amplification of small inhomogeneities. Statistical characterization of clustering is important 
for understanding the dynamics of gravitational clustering. Correlation functions and count in cell statistics 
are among the oldest and most widely used statistics to quantify clustering. Evaluation of higher order 
correlation functions with full shape dependence is difficult given limited size of present day galaxy catalogs 
and numerical simulations (Fry & Peebles 1978, Peebles 1980). Much progress has therefore been made 
in estimation of volume averages of these quantities (Gaztanaga 1992, Bouchet et al. 1993, Gaztanaga & 
Frieman 1994, Colombi et al. 1996a, Szapudi et al. 1996), which can be directly related to moments of one 
point count probability distribution function. Analytical methods based on tree level perturbation theory and 
loop corrections were used to evaluate these quantities in the quasi-linear regime (Peebles 1980, Juszkiewicz 
et al. 1993, Bernardeau 1992, Bernardeau 1994, Munshi et al. 1994, Colombi et al. 1992, Scoccimarro & 
Frieman 1996a-b,). In the highly nonlinear regime, although no full theory exists, scaling models are often 
used to predict their behaviour (Peebles 1980, Balian & Schaeffer 1989). Computational methods were also 
developed based on factorial moments to subtract Poisson noise from discrete data. Errors associated with 
these one point cumulants and other related quantities such as void probability function have also been 
estimated and correction procedure have been developed (Colombi et al. 1995, Szapudi & Colombi 1996, 
Munshi et al. 1997). 

Although one point quantities carry useful information about dynamics and background geometry av- 
eraging associated with such quantities causes a significant loss of information. Alternative methods use 
correlation of pairs of cells (Szapudi et al. 1992, Meiksin et al. 1992, Szapudi et al. 1995). Cumulant corre- 
lators provide a natural generalization of one point cumulants. The hierarchal ansatz in the strong clustering 
regime and perturbation theory in the weak clustering regime have definite predictions for these quantities. 
Using cumulant correlators it is possible to separate amplitudes associated with different tree topologies up 
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to fifth order in the correlation hierarchy. Szapudi & Szalay (1997) have analysed APM data using cumulant 
correlators. We use these statistics to study scale-free simulations in two and three dimensions (2D and 3D). 

In the next section we outline the theoretical framework necessary for using cumulant correlators, pre- 
dictions from perturbation theory, hierarchal ansatz and extended perturbation theory. In section §3 we 
describe our numerical simulations and techniques to evaluate cumulant correlators from simulation data. 
We discuss implications of our results in section §4 and compare it with existing results from galaxy catalogs. 



2. Theoretical Predictions 

We define factorial moment correlators of a pair of cells of volume I 3 separated by a distance r as 



((jVi) fc (^ 2 )z)-((jVi) fc )((jV2)z; 
<iV>* 



w k i{l,r) = /AAfc+i > k^QJ^O, (1) 



and the normalized one point factorial moment 

w feo(0 = {Ni)k , (2) 

where we have used the notation (N)k = N(N — 1) . . . (N — k + 1), and the angular braces (. . .) denote 
average over all possible positions of the cells. 

It is possible to relate Wko and wu to normalized one point cumulants Sn = {S N ) / (<5 2 )^ -1 ^ and cumulant 
correlators Cnm = (S N \xi)S M (xa))/ (5(xi)5(x2)}(5 2 } < ' N+M ~ 2 \ However we find it to easier to work with 
Q N = S N /N<> N -^ a,ndQ NM = C N m/N^ n ~ v >M^ m ~ x ^. While central moments can also be used to computed 
these quantities, factorial moments are known to be more suitable for subtracting Poisson shot noise. 

One can introduce the generating function for the factorial moments in terms of the cumulants Qn- 

oo 

W(x) = exp r N x N Q N (3) 

where we have defined T N = N N - 2 £ s N ~ l /N\ and £ s is the volume average of £2 over volume of the cell i 3 . 

The generating function can be linked with the one point void probability distribution function. In 
general Qn parameters show scale dependence, increasing from quasi-linear phase to highly nonlinear phase. 
Hierarchal form of correlation functions demand Qn to be independent of scale in the nonlinear regime. 

We can also construct generating function of factorial moments W(x,y) = W mn x M y N / mini which 
can be related to generating function of cumulant correlators Q(x, y) 

00 

Q(x,y)=£ l x M y N Q NM T M T M NM, (4) 

M=l,iV=l 

by the following equation 

W(x,y)=W(x)W(y)(e X pQ(x,y)-l). (5) 
Expanding equation (ph one can compute cumulant correlators Qnm from factorial moments Wnm- 



Q\2^l^2£,r - 

)i33rir 3 £ r = W13/6 - 

Q 22 r 2 ; 4£ r = W22 /4 



W12/2 - £r 
W12/2 - W20/2 + Cr 
- W12 + tr~ C 2 /2 



(6) 
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As it is clear from equations (^J), cumulant correlators are determined by two different length scales, 
length scale corresponding to £ s and scale corresponding to £ r . We use the terms quasi-linear and nonlinear 
depending upon values of £ s . 

2.1 Nonlinear Regime 

The extreme nonlinear stage of gravitational clustering is believed to be well described by a power law 
solution for correlation function £2 = T ~ 7 , where 7 can be expressed as a function of initial power spectrum 
(Davis & Peebles 1977, Balian & Schaeffer 1989. All higher order correlations exhibit a self similar scaling 
in strong clustering regime 

6v(Arx, . . . Ai-at) = A-^-^Mri, ■ ■ ■ r N ). (7) 

More explicit functional forms for TV-point correlations are constructed by summing over products of N— 1 
two-point correlation functions corresponding to different topologies, each of which representing a tree graph 
spanning N-points with amplitudes T n . a 

a, N— trees labelings edges 

Tree graphs spanning three points can have only one topology. The amplitude associated with this graph 
Q3 contributes through three different configurations of these points, so that cumulant correlators to have 
following form: 

(<5(x 1 ) 2 5(x 2 )) = 2Q 12 Ur = Q 3 (2Ur + &). (9) 

Higher order cumulants get contributions from different types of tree diagrams. At fourth order trees 
connecting four points have two different topologies, known as the snake (with its amplitude denoted by R a ) 
and the star topology (corresponding amplitude denoted conventionally by Rb). Using this notation we can 
write 



<<$l(xi) 3 <J 2 (x 2 )> = 9Qi 3 &.£ = ^r^Ra + SHrgRb + 6£,%Ra + &Rb (10) 
( ( 5(x 1 ) 2 5(x 2 ) 2 ) = 4Q 22 ^ 2 = AUlRa + ^rtsRb + ^AsRa + 4&Rt. (11) 

These simultaneous equations are linear in R a and Rb and can be solved once Q22 and Q13 are determined. 
When £ r < < 1 which will be the case when two cells are far away, one can define linear cumulant correlators 
(linear in £ r ) by considering only terms linear in £ r . The linear solution to equation (|ll|) are R a = Q22 and 
Rb = 3Qi3 — 2Q 22 . In such a linear accuracy one can show that 

Qnm ~ Qn+m (12) 

Although amplitudes of tree terms with different topologies can be estimated by cumulant correlators, the 
situation becomes more complicated at higher order as number of degenerate correlators become less than 
number of topologies, making the system indeterministic. Other interesting questions regarding scaling of 
generating functionss Q{x,y) and related question about the nature of Qnm in the highly nonlinear regime 
can be addressed when bigger N-body simulations become available. 

2.2 Quasi-linear Regime 

In quasi-linear regime when £ s is smaller than unity it is possible to expand 5 = S^ 1 ' + 5^' + S 1 - 3 ^ . . . where 
perturbation expansion is valid as long as the series is convergent. Using perturbative calculation of two 
point quantities it was possible for Bernardeau (1996) to express C pq at lowest order. Perturbative terms 
contributing to the expansion of (5 p (xi)<5 9 (x 2 )) can be written as 
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(^( Xl )^(x 2 )) « J2 (n^ i) (xi)f[^ ) (x2)> (13) 

decompositions i— 1 J — 1 

where sum is taken over all possible decompositions and pi and qi satisfies 

p <? 

E^ + E*=2(P + 9)- 2 - (14) 

i=l .7=1 

In lowest order in perturbation theory it is possible to write (Bernardeau 1996), 

<<5 p (x 1 )5'(x 2 )) c = ^(^(x))^- 2 ^!)^)) = C pq ^-^ r . (15) 

Using method of generating function it was shown that C pq can be decomposed into following relation 
(Bernardeau, 1996), 

C pq = CplCql. (16) 

In 3D the lowest order C pq can be expressed in terms of spectral index n, 

C 2 i =68/21-l/3(n + 3) (17) 
C31 = 11710/441 - 61/7(n + 3) + 2/3(n + 3) 2 . (18) 

It should be noted that such factorization is possible only in case of tree-level diagrams also calculations 
were done using large separation limit. Given the whole hierarchy of Cnm it is possible to compute bias 
associated with gravitational clustering in the quasi-linear regime (Bernardeau 1996). 

2.3 Connecting different regimes: Extended Perturbation Theory 

Tree level perturbation theory can predict Sn parameters for top-hat smoothing, it was used (as out lined 
above) to compute cumulant correlators Qn,m to arbitrary order in large separation limit. The expressions 
derived for cumulants depends on index of initial power spectrum. However it was later realized that the same 
quasi-linear expressions can describe evolution of Sn parameters from quasi-linear regime to the nonlinear 
regime by allowing spectral index n to vary with scale (Colombi et al. 1996b). It is interesting that whole 
hierarchy of Sn can be described by such a variation of effective spectral index n e f / (unfortunately no straight 
forward relation exists between true nonlinear power spectra and n e ff). An extension of perturbation theory 
to the non-perturbative regime was considered by Szapudi & Szalay (1997) for cumulant correlators Qn,m- 
Tree level perturbation theory as mentioned before predicts C pq — C v \C q \. Using same logic as extended 
perturbation theory one can expect such a relation to hold even for small cell sizes at large distances although 
separately each of these term may vary significantly. 

3. Simulations and Data Analysis 

The simulations used here are numerical models for the gravitational clustering of collisionless particles in an 
expanding background. We study evolution of initial Gaussian perturbations in £1 = 1 universe. All the 2D 
simulations are done with a particle-mesh (PM) code with 512 2 particles with an equal number of grid points 
and in 3D 128 3 particles with 128 3 grid points The code has at least twice the dynamical resolution of any 
other PM code with which it has been compared. The 2D simulations are described in detail in Beacom et. 
al (1991), with a video of the evolution in Kauffmann and Melott (1992). The 3D simulations are described 
in Melott and Shandarin (1993). Both sets, which consist of multiple realizations (different random seeds) 
for a variety of power spectra, have been widely used for comparative studies of various statistical methods 
and dynamical approximations. 
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Fig. 1 . — Lower panel displays factorial moment correlators Wki as a function of separation r in 2D in units 
of grid scale. Degenerate parallel lines correspond to increasing values of TV + M from 2 to 6. Left panels 
display results for scale free spectra with power law index n = — 2, middle panel for n — and right panel 
n = 2 in two dimension. Middle and upper panels plot cumulant correlators of different orders. Dashed 
curves in middle panels represent QniQmi for different values of N and M, closest solid curves represent 
Qnm for same values of N and M. According to predictions from extended perturbation theory they should 
overlap Qnm = QniQmi- Top most panels compare predictions of hierarchal ansatz Qnm — Qn+m- 
Adjacent curves correspond to different N and M with same N + M. Error bars were calculated by finding 
scatter in results with different realization of same initial power spectra. 
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Fig. 2. — Same as Figure - 1 but results from only one realization are plotted to show level of agreement 
with theoretical predictions in individual realizations. 
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In this paper we analyze a subset of the simulations with featureless power-law initial spectra of the 
general form, 

P(k) oc k n for k < kc, (19) 
= for k > k c . (20) 

We have analyzed power-law models with n = 2, 0, —2 in 2D and n = 1, 0, —1, —2 in 3D with a cutoff in each 
case at the Nyquist wave number: k c — 256 kf for 2D and k c — 64 kf for 3D where kf = 27r/ibox is the 
fundamental mode associated with the box size. 

We choose ct(/cnl), the epoch when the scale 2-7t/A:nl is going nonlinear as a measure of time. 

-1 



; fc/ NY ^(fc)fck 
J^P(k)kk^ 



The first scale to go nonlinear is the one corresponding to the Nyquist wave number. This happens, by 
definition, when the variance a is unity. Of course as a increases successive larger scales enter in the nonlinear 
regime. The simulations were stopped at Anl = 2^ gr ;d, 4Z gr id, 8Z gr id, Lb ox /2. In our study in 2D we took 
the epoch when L^ ox /lQ for analysis of n = — 2 spectrum and for n = and2 we studied the epoch when 
Lbox/8 is going nonlinear. In 3D we had less dynamic range so choose the epoch when Lbox/4 for analysis 
of n = — 2 spectrum and for n — —1,0, andl we studied the epoch when i box /2 is going nonlinear. The 
separation r that we study is much less compared to the length scale going nonlinear so it is unlikely that 
our results will be affected much by boundary conditions. 

The growth rate of various modes in the linear regime were studied by Melott et al. (1988) for this PM 
code. The results at A — 3Z gr id are equivalent to the ones obtained by a typical PM code at A = 8Z gr id, due 
to the staggered mesh scheme. So we expect that our code performs well at the wavelength associated with 
four cells and since the collapse of 4 l gr id-size perturbations will give rise to condensations of diameter 2 Z gr id 
or less, the smallest cell size that can be safely resolved is 2Z gr id- 

Computations of count probability distribution function (CPDF) Pi(n) and two-point count probability 
distribution P; jT -(n, m) were done by laying down a grid of mesh spacing I and counting occupation number in 
each cell to compute the probability of finding n objects in cell size I and also the joint probability of finding 
n and m objects in two cells separated by a distance r. Statistics were improved by perturbing the grid in 
each orthogonal direction and keeping the mesh undistorted while repeating the counting process. Both in 
two and three dimensions we considered only cells of size l gr id- While this particular choice of cell size is open 
to question, we will show that our results match remarkably well with our earlier studies done using bigger 
cell sizes. With this cell size we could reach probabilities as few times 10~ 6 in 2D and 3D. Computations of 
factorial moments and factorial moment correlators were done using computed values of Pi{n) and Pi, r {n, m). 
Computed values of £2 for cells used were found to be 8.98, 27.87, 60.39 in 2D for spectral index —2, 0, and 2 
respectively. In 3D these values were found to be 25.46, 88, 125.5, and 171 respectively. 

Different spurious effects such as shot noise and finite volume corrections are important while computing 
count probability distribution functions. For small cell sizes Poisson noise starts dominating, whereas larger 
cell sizes are dominated by finite volume corrections. We find Wki parameters to be dominated by shot noise 
with increasing separation of cells, this effects starts to be severe with higher order cumulants. This effect 
starts dominating even for small separation for spectra with less power on larger scales. Using large cells 
produces effect similar to smoothing the density field and hence reducing the effective correlation length 
scale between them. These two dominant effects reduce the range of separation which can be probed for 
studying cumulant correlators. Finite volume corrections are more difficult to quantify. Methods based on 
scaling of count in cell statistics were shown to be effective in correcting finite volume effects (Munshi et al. 
1997), similar arguments can be used to develop corrections for cumulant correlators. The validity of such a 
method depends on correctness of the scaling model, which we test in this paper. Developments of methods 
to correct Qnm for finite volume effect are left for future work. 

Computed Wki for 2D simulations are plotted in Figure - 1 as a function of cell separation r. Corresponding 
results for 3D are presented in Figure - 5. Parallel degenerate lines correspond to different values of TV + 
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Fig. 3. Lower order hierarchal amplitudes Q3, R a , Rb calculated from the fully nonlinear cumulative 
correlators in 2D arc plotted against separation r measured in grid units. Solid curve represent estimates 
of Q3 lower and upper dashed lines represent amplitudes of fourth order snake and star graphs R a and Rb 
respectively. Open triangles and square correspond measurement of Q3 and Q4 respectively from factorial 
moments of counts in cell statistics after finite volume corrections were taken into account (from Munshi 
ct al. 1997). Filled triangles and squares are measurements of same quantities using cumulant correlators 
without finite volume correction. 



Fig. 4. — Each panel displays factorial moment correlators Wki as a function of separation r in unit of grid 
scale for different initial scale free power law spectra in 3D. Scatter in each plot is computed from scatter 
in four different realizations of same power spectra. 
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Fig. 5. — Nonlinear cumulant correlators is plotted against separation r in 3D as measured in units of grid 
scale. Uppermost curve correspond to Q31, lowest curve to Q22 and the dashed one to Q\ x - Solid and dashed 
straight lines correspond to predictions for Q31, and from perturbations theory at large separation. 



Fig. 6. — Hierarchal amplitudes Q3, R a , Rb calculated from the nonlinear cumulative correlators in 3D 
are plotted against separation r measured in grid units. Solid curve represent estimates of Q3 lower and 
upper dashed lines represent amplitudes of fourth order snake and star graphs R a and Rb respectively. Open 
triangles and square correspond to measurement of Q3 and Q4 respectively from factorial moments of counts 
in cell statistics after finite volume corrections were taken into account (from Munshi ct al. 1997). Filled 
triangles and squares are measurements of same quantities using cumulant correlators without finite volume 
correction. 
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M increasing from bottom to top in each panel. Middle and top panels show variation of Qnm with r. 
Middle panels test the validity of extended perturbation theory. Dashed lines correspond to QniQmi while 
neighboring solid lines represent Qnm for same values of N and M. Error bars are computed by estimating 
the scatter in results from four different realization of each spectra. It is interesting to note that while Qnm 
shows an increasing trend with increasing N + M for spectra with large scale power n = —2 in 2D, the 
trend is reversed for n = and 2. The results presented in 2D are for N + M = 4, 5 and 6 respectively. In 
3D meaningful computation of Qn,m were possible only for N + M < 4. We find that computed values of 
Q nm are more stable towards fluctuation at large separation if one of the indices of cumulant correlators is 
larger than the other (e.g. Wsi is more stable at larger separation as compared to 1022)- Close associations 
of dashed curves with solid curves proves the validity of extended perturbation theory. It is also easy to 
notice that agreement is better in case of 2D compared to 3D. Larger simulation size might be a possible 
explanation for such an effect. Given that we restrict our result to highly nonlinear regime, it is intcrsting 
to note that computed values of Cnm are not completely different from predictions of perturbation theory 
at large separation. In Figure - 5, dashed and solid straight lines are predicted values of Q21 and Q31 from 
equation (17) and equation (18). Top panels in Figure - 1 and Figure - 5 tests the validity of equation (|l2| ) in 
2D and 3D respectively. Nearest neighbor curves correspond to different values of N and M which produces 
same N + M. Agreement in this case is similar to that of extended perturbation theory. 

Using cumulant correlators it is possible to compute amplitudes of different tree topologies contributing 
to four and five point correlation functions. We have computed amplitudes Q, R a and Rb- Figure - 3 shows 
that these amplitudes are almost constant in 2D in the highly nonlinear regime. In 3D Figure - 6 shows large 
fluctuations in computed values of Rb although in general they are fairly constant within the limited range of 
nonlinearity studied by us. Open triangles and squares denote measured values of Q3 and Q4 using factorial 
moments of CPDF. On the other hand solid triangles and squares denotes measured values of Q3 and Q4 
from cumulant correlators. Q4 was calculated by using the relation I6Q4 = 12i? a + 4i?j,, from computed 
average values of R a and Rb. The remarkable agreement of both the methods increases confidence in results 
based on cumulant correlators. It is also to be noted that whereas computation of one point cumulants were 
done by taking volume corrections into account, such effects were neglected in the computation of cumulant 
correlators. Our results show that Q4 is much closer to Rb which might indicate that star topologies in general 
dominate over snake topologies. However, more systematic studies are necessary which will incorporate finite 
volume effects and have more dynamic range than this study. 

6. Discussion 

Measurements of cumulant correlators and lower order one point cumulants have already been done for Lick 
catalog, SDES and the APM survey. Analyzing cumulant correlators in the APM survey, Szapudi & Szalay 
(1997) demonstrated the validity of the hierarchal ansatz to unprecedented accuracy. They found Q3 = 1.15, 
R a = 5.3 and Rb = 1.15. These results were also shown to be in good agreement with computation of 
one point cumulants based on factorial moments Q3 = 1.7 and Q4 = 4.17. Measurements form SDES give 
Q3 = 1.16 and Q4 = 2.96. It was suggested that slightly low value of Q4 form SDES was lack of nonlinear 
corrections in their computations. Other measurements from APM produces Q3 = 1.6 Szapudi et al. (1996) 
and Q3 = 1.7 by Gaztanaga (1994). The values of fourth order cumulants found by these authors were 
Q4 = 3.2 and Q4 = 3.7 respectively. 

Use of methods based on cumulant correlators by Szapudi and Szalay has produced lower values for Rb 
compared to R a . In our simulations both in 2D and 3D we find that relative position of these two amplitudes 
depend on initial power spectra. In all cases Rb was found to be either equal to or greater than R a , although 
the separation was found to decrease with increasing n. For power spectra such as n = 1 in 3D and n = 2 
in 2D which have more power at smaller scales these two quantities arc found to almost coincide with each 
other. Although this is in disagreement with findings of Szapudi and Szalay (1997), our results seems to 
be closer to that of by Fry & Peebles (1978) who found R a = 2.5 ± 0.6 and R b = 4.3 ± 1.2 from their 
analysis of Lick catalog. At any rate, they are analyzing data in projection, and we a simulation with full 
information. We do not know to what extent projection and the possible inclusion of non-gravitational 
physics may influence the results. 
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There have been several attempts to solve BBGKY equations in the highly nonlinear regime, although 
no general solution exists so far. Attempts were made to solve these equations by assuming specific sep- 
arable hierarchal forms for higher order correlations in phase space. Fry (1982) derived Qn,u = Qn — 
{4Q/N) N ~ 2 N/(N - 2), in particular this model predicts R a = R b = Q 4 = 2Q 2 /3. Hamilton (1988) redoing 
the analysis found only snake topologies to contribute in higher order correlations while all non-snake con- 
tributions vanish identically which gives Qn — (Q /N) N ~ 2 Nl/2. For lower order amplitudes this predicts 
Ra = Q 2 1 Rb = and Q4 = 3Q 2 /4. The ansatz used by Bernardeau and Schaeffer (1992) predicts R a — Q\. 
While all these efforts to solve BBGKY equations provide useful insight into complications arising due to 
the nonlinear nature of the problem, it is clear from our study that cumulant correlators are a very useful 
tool in comparing such predictions with computer simulations. The general agreement of 2D and 3D results 
reinforces the many other results which have shown that 2D may be a useful tool to extend the dynamic 
range for exploration of statistics in cosmology. 

Acknowledgements: D.M. acknowledges support from PPARC under the QMW Astronomy Rolling Grant 
GR/K94133. A.L.M. wishes to acknowledge the National Center for Supercomputing Applications for sup- 
port to perform the ensemble of simulations, and the financial support of the NSF EPSCoR program. 
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